function I = romberg(f,n)
m=floor(log2(n))+1;
T=zeros(m,m);
for i = 1:m,
    T(i,1)=trapezoid(f,2/(2^(i-1)));
    for j=2:i
        T(i,j)=((4^j)*T(i,j-1)-T(i-1,j-1))/((4^j)-1);
    end
end
I=T(m,m);
